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ABSTRACT 

We report results of 3D MHD simulations of the dynamics of buoyant bubbles in magnetized galaxy 
cluster media. The simulations are three dimensional extensions of two dimensional calculations 



reported by Jones & Dc Young ( 2005 ) . Initially spherical bubbles and briefly inflated spherical bubbles 
all with radii a few times smaller than the intracluster medium (ICM) scale height were followed as 
they rose through several ICM scale heights. Such bubbles quickly evolve into a toroidal form that, 
in the absence of magnetic influences, is stable against fragmentation in our simulations. This ring 
formation results from (commonly used) initial conditions that cause ICM material below the bubbles 
to drive upwards through the bubble, creating a vortex ring; that is, hydrostatic bubbles develop into 
"smoke rings", if they are initially not very much smaller or very much larger than the ICM scale 
height. 

Even modest ICM magnetic fields with (3 = P gas / 'P mag < 10 3 can influence the dynamics of the 
bubbles, provided the fields are not tangled on scales comparable to or smaller than the size of the 
bubbles. Quasi-uniform, horizontal fields with initial (3 ~ 10 2 bifurcated our bubbles before they rose 
more than about a scale height of the ICM, and substantially weaker fields produced clear distortions. 
These behaviors resulted from stretching and amplification of ICM fields trapped in irregularities 
along the top surface of the young bubbles. On the other hand, tangled magnetic fields with similar, 
modest strengths are generally less easily amplified by the bubble motions and are thus less influential 
in bubble evolution. Inclusion of a comparably strong, tangled magnetic field inside the initial bubbles 
had little effect on our bubble evolution, since those fields were quickly diminished through expansion 
of the bubble and reconnection of the initial field. 

Subject headings: galaxies: clusters: general - MHD- galaxies: active galaxies 



1. INTRODUCTION 

There is now ample evidence from radio and X-ray 
observations that active galactic nuclei (AGN) gener- 
ate energetic structures that continue to interact with 
galaxy cluster environments even after the central engine 
shuts down. Observations of detached bubbles of radio 
plasma and coincident X-ray cavities in cluster cores il- 
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mated from X-ray observations, which suggest that up- 
wards of 10 59 — 10 60 e rgs of energy are present i n these 
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The evolution of relic bubbles is an interesting problem 
in part because their structures are seen to remain coher- 
ent longer than analytic estimates would suggest. The is- 
sue of bubble fragmentation and its associated timescale 
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is of importance for models that use "AGN feedback" 
to reheat the ICM and thus suppress both ICM cool- 
ing flows and the occurrence of large amounts of star 
formation in the cluster core. This is because such re- 
heating needs to be spread throughout the ICM in the 
central regions, and if radio bubbles do not fragment on 
timcscalcs comparable to the local ICM cooling time, the 
distribution of their energy over large volumes requires 
the assumption of some other, yet unspecified mecha- 
nism. This fragmentation question is also related to the 
more general problem of AGN feedback in current cos- 
mological models and in particular to the use of "radio 
AGN feedback" to produce the observed distrib ution of 
massive ga laxy morphologies and colors (e.g., |Croton| 
et al.|2006|. In particular, because the bubbles are light, 



presumably being filled with very hot and possibly rela- 
tivistic plasma, they buoyantly rise in the cluster poten- 
tial and are subject to disruption by both the Raylcigh- 
Taylor (R-T) and Kelvin-Helmholtz (K-H) instabilities. 
Simple linear stability analysis suggests that disr uption 
should take place o n timescales of ~ 10 7 years (see Heinz 



& Churazov 2005 for example), yet intact bubbles are 
seen at distances from the cluster core t hat imply rise 
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times of an order of m agnit ude longer ( Birzan et al 
H2005[ ). As |De Young] ( |2003p pointed 
ICM magnetic fields and possibly in- 



possibly 

ternal bubble magnetic fields could play an important 
role in relic evolution by stabilizing bubble structures 
against these instabilities for timescales ~ 10 8 years. On 
the other hand, the role of a magnetic field in either of 
these instabilities is dependent on both the strength and 
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the orientation of the field with respect to the unstable 
boundary. In addition, the field does not enhance sta- 
bility for modes orthogonal to the plane containing the 



field a nd the boundary (e.g., Chandrasekhar[l961 Dursi 
20071. Furthermore, nonlinear evolution of R-T and K 



H instabilities can be very difficult to predict, even in 
the hydrodynamic case, but especially in the presence of 
magnetic fields; disruptive trends can be both enhanced 
and reduced compar ed to linear pred i ctions, dependin g 



on the details (e.g., Jun et al. 1995 Ryu et al. 2000 



Consequently, understanding of this issue requires fully 
nonlinear MHD numerical simulations. 

Several groups have previously conducted simulations 
in 2D and 3D to explore basic bubble dynamics and mor- 
phology and t heir relationship t o ICM magnetic fields. 



For example, Churazov et al. 
buoyant rise times of ~ 10 s 



(|2001| first established 
years for bubbles in clus- 
ter environments in 2D hydrodynamic simulations. Hy- 



drodynamic simu lations conducted by|Briiggen fe Kaiser| 
( |2001| in 2D and |Briiggen et al.| ( [2002P in 3D confirmed 
the anticipated timescales of ~ 10' years for disruption 
from hydrodynamical instabilities in such bubbles. To 
explore the effec ts of magnetic fields o n bu bble dynam- 



ics and stabil ity, |Robinson et al. ( 2004 1 and Jones & De 
Young ( |2005[ ) conducted multiple series of 2D magne- 
tohydrodynamic (MHD) simulations involving relatively 
homogeneous ambient magnetic fields with a variety of 
field orientations with respect to the cluster gravitational 
field. Once again the bubble rise times were ~ 10 s years. 
These calculations also demonstrated that even ambient 
fields initially much too weak to provide linear stability 
to R-T and K-H instabilities could, nonetheless, appar- 
ently prevent bubble disruption on timescales compara- 
ble to the inferred lifetimes of the observed "relic" ra- 
dio bubbles in clusters. The surprising ability of such 
weak fields to stabilize bubble boundaries in those simu- 
lations came from the stretching of field lines in response 
to the bubbles' rise through the ICM and from vorti- 
cal motions generated within th e bubbles. More recent 
work in the 3D MHD regime bylRuszkowski et al.|(|2007|) 
has suggested, on the other hand, that random ambient 
magnetic fields only provide such bubble stability if the 
coherence lengths of these fields are large with respect to 
the bubble size. 

Magnetic fields may also play significant roles in the 
dynamics and thermal conduction properties of "cold 
fronts", which are contac t discontinuities separating 
merging subclusters (e.g. 



Vikhlinin et al. 
et aL||2007| iDursi & Pfrommer||2008p. That 
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physics problem in some respects, although strong dif- 
ferences in the initial conditions, relative properties of 
the media, and the overall dynamics make direct com- 
parisons between magnetic field roles in cold fronts and 
bubble boundaries nontrivial. We will not address mag- 
netic effects in cold fronts here. 

In this work, we explore further the three-dimensional 
evolution of bubble dynamics and morphology in the 
presence of magnetic fields. We present results from a se- 
ries of 3D MHD simulations of buoyant bubbles i nflated 
in a cluster atmosp here similar to that used by |Jones 
fc De Young| ((2005). We test, for several field strengths 



and topologies, what determines bubble shapes and how 
magnetic fields affect the evolution of bubbles. Since 
the magnetic fields in the bubbles have distinct origins 



from the cluster fields they encounter, we take pains in 
our initial conditions to isolate the bubble and ambient 
medium magnetic fields. We also explore the dynami- 
cal evolution of buoyant bubbles that are inflated by the 
injection of energetic material into an initially small vol- 
ume. Many of the previously published calculations were 
based on bubbles that were initialized as static spher- 
ical structures inside an ICM. Since real radio bubbles 
are remnants of AGN outflows and thus not likely to be 
particularly spherical at the start, it is important to un- 
derstand how simulation results depend on this choice of 
initial conditions. 

The paper is organized as follows. Section [2] describes 
the details of the calculations performed, including a de- 
scription of the numerical algorithm and explanations of 
the simulated models. Section [3] contains a discussion of 
the results and an analysis of the physics of bubble prop- 
agation and magnetic field evolution. Conclusions and a 
summary are provided in § |4j 

2. CALCULATION DETAILS 

The compu tational approach we desc ribe here is simi- 
lar to that of | Jones fc De Young (2005) but extended to 
three dimensions. We inflated bubbles from a spherical 
volume of radius rj, near the base of a plane stratified 
ICM atmosphere initially in hydrostatic equilibrium. 

In most of the simulations the ICM was magnetized 
with magnetic field pressures that were generally much 
smaller than the ICM thermal pressure. In one set of sim- 
ulations, the ambient magnetic field was quasi-uniform, 
meaning that the field asymptotically achieves uniform 
orientation far from the bubble inflation region. In an- 
other set of models, the ICM field was tangled on scales 
either comparable to or larger than the bubble inflation 
region. Cases are included with magnetic field strengths 
interior to the inflated bubbles comparable to the exter- 
nal fields and with bubble fields absent. For comparison, 
we also include two strictly hydrodynamic cases without 
any ICM magnetic fields. General simulation descrip- 
tions follow, with specific details of t he ICM given in § 
PET] and those of the bubbles in § [2~2| 

Our simulations were conducted on a 3D Cartesian 
grid of total physical dimensions x S i ze = 40 kpc, y S ize = 
z s ize = 34 kpc, covering the ranges x = [5 kpc : 45 kpc], 
y = z w [— 17 kpc : 17 kpc]. The gravitational accel- 
eration is in the —x direction in this coordinate sys- 
tem. The bubble origin was centered at Xb = 13 kpc, 
yi, — Zb — 0. The grid was partitioned into uniform com- 
putational zones of size Ax = Ay — Az — 133 pc. At 
the y and z boundaries, continuous boundary conditions 
were maintained. To maintain hydrostatic equilibrium 
in the x direction, extrapolation boundaries were em- 
ployed at x = 5 and x = 45 kpc. In particular, the 
density in a given boundary zone was set through a lin- 
ear extrapolation of the densities in the adjacent pair of 
physical grid zones. Using this updated boundary den- 
sity, the boundary pressure was derived assuming hydro- 
static equilibrium within the boundary zones. Addition- 
ally, the velocities were chosen to be continuous across 
both x boundaries, in order to minimize the impact of 
waves incident upon those zones. The calculations were 
terminated when any bubble material approached any 
boundary of the computational domain (generally x top ). 

Our simulations employed a second-order, nonrelativis- 
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TABLE 1 

Summary of Simulations of Buoyant Bubbles 



ID 1 ICM Field Geometry a 2 B a (fiG) 23 U (Myr) t e (Myr) 



US 


Uniform 


120 


5 





150 


US-I 


Uniform 


120 


5 


10 


136 


US-F 


Uniform 


120 


5 





150 


UM 


Uniform 


3000 


1 





150 


UM-I 


Uniform 


3000 


1 


10 


111 


uw 


Uniform 


oo 








150 


UW-I 


Uniform 


oo 





10 


132 


TlS 


Tangled-L 


120 


5 





108 


T L S-I 


Tangled-L 


120 


5 


10 


111 


T L S-F 


Tangled-L 


120 


5 





108 


T B S 


Tangled-B 


120 


5 





114 


T S S-I 


Tangled-B 


120 


5 


10 


130 



* The first letter indicates the ambient magnetic field geometry: uniform (U) or tangled (T) with 
coherence lengths that arc long (sub L) or bubble-sized (sub B). The second letter refers to the strength 
of the field: Strong (S), Moderate (M), or Weak (W). Letters to the right of the dash indicate cither that 
the bubble is inflated (I) or that the bubble itself contains a magnetic field (F), while models without 



these indicators imply the opposite case (i.e., non-inflated or without a bubble field), 
true in a statistical sense for the tangled fields. 3 Measured at x — 30 kpc. 



This is only 



tic, Eulerian, total variati on diminishing (TV D), ideal 
3D MHD co de, described in |Ryu fc Jones ( 1995 1 and Ryu 
et al. ( 1998 1. The code explicitly enforces the divergence- 
free condition for magnetic f ields through a co nstrained 
transport scheme detailed in Ryu et al. (19981. The nu- 



merical method conserves mass, momentum, and energy 
to machine accuracy. To set up hydrostatic equilibrium 
in the ambient medium in the simulations presented here, 
acceleration due to the assumed gravity was added as a 
source term in the —x direction through operator split- 
ting. The corrected x momentum at each timestep was 
used to recalculate the total kinetic plus thermal and 
magnetic energy. In this simplified treatment of gravity, 
momentum and energy are no longer exactly conserved, 
but we have verified that associated errors are much too 
small to influence our results. A gamma-law gas equa- 
tion of state was assumed with 7 = 5/3. A passive mass 
fraction or "color" tracer, C/, was introduced for ma- 
terial originating inside the bubble inflation region. C/ 
was set to unity inside the bubble, while C j = in the 
ambient medium. 

2.1. Gravity and The Ambient Medium 
We applied the same gravity model as |Jones fc De 



Young ( 2005 1 , which was derived from the combined 
mass distributions of the central active galaxy and the 
cluster. The cluster had an NFW mass profile, while the 
active galaxy had a King model. The galaxy and cluster 
mass normalizations were chosen such that the cluster 
mass inside 10 kpc was 3.5 x 10 10 M Q and the galaxy 
mass inside 20 kpc was 3.5 x 10 12 M Q . Inside r = 60 kpc 
(i.e., everywhere on our computational domain) matter 
associated with the active galaxy dominated the gravita- 
tional force. 



To match t he planar ICM geometry of Jones & De 
Young (20051, we projected the resulting gravitational 
acceleration onto the x direction; that is, the gravity 
was one dimensional, leading to a plane stratified hy- 
drostatic equilibrium for the initial ICM. That ICM was 
assumed to be isothermal with T 3.5 x 10 7 K, or, 
cquivalently, kTo = 3 keV. We assumed a hydrogen mass 
fraction X = 0.75 (for a mean molecular weight /i = 0.6), 



giving an adiabatic ICM sound speed c Si o = \/jPoJpo = 
894 km s -1 = 0.9 kpc Myr -1 . The electron density was 
set to n e = 0.1 cm~ 3 at x = 5 kpc, leading to the atmo- 
spheric density, p , and pressure, Pq, profiles shown in 
Figure]!] Th ey are identical to those in the IJones & De 
Young ( 2005 1 simulations. The scale height of the ICM, 
here defined to be h — c 2 /\g\, turns out to behave as 
h « x, which allows us to approximate the density and 
pressure distributions as a simple power law for analysis 
purposes; in particular 



Po(x) 
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Po(x) = P x 



(1) 



where 7 = 5/3 is the gas adiabatic index, xq is some 
fiducial vertical position, and P Xo is the pressure at xq. 

We explored consequences of two simple ICM mag- 
netic field models, one quasi-uniform and one tangled, 
with both designed so that the ICM field did not ini- 
tially penetrate into the bubble. The quasi-uniform ICM 
field, which we identify using the letter "U" in model 
labels, was a 3D extension of the 2D field model used 



by Jones & De Young ( 2005[ ). In this case field lines 
initially lay in the x — y plane, being derived from a 



vector potential A z (r, 



B a r[l - (r b (z) /r) 2 ] cos </>, de- 



fined in cylindrical coordinates (r, (f>, z) with respect to 
the center of the bubble inflation region with </> = along 
the x direction (i.e., the direction of bubble propaga- 
tion). In this expression r\,{z) — — (Jz — Zb)/ r b) 2 
for \z — r&| < rfc. At large distances from the bubble, 
the field lines in this model become horizontal, but they 
are tangent to the bubble surface when r = rb(z), so 
that they skirt over its top and bottom in the vertical 
plane. Figure [2] shows the geometry of the quasi-uniform 
ICM field. The field strength parameter Bo was scaled 
vertically as B (x) oc yJP (x), so that the ratio of gas 
to magnetic pressures, /3 = 8ttP /Bq was constant with 
height away from the inflation region. Thus, the ambient 
magnetic field strength generally decreased with height 
in the ICM. 

To provide a direct co mparison to the 2D models of 
Jones & De Young ( 2005 1 we carried out in the present 




Fig. 1. — Initial buoyant bubble environments. Top left: Gravitational mass as a function of radius from the core: total (solid line), galaxy 
(dotted line), and cluster (dashed line). Top right: Gravitational acceleration. Bottom left: Gas density. Bottom right: Gas pressure. 



study several quasi-uniform, "U" , field model simulations 
with the same values of (3q that they chose; namely, /3q — 
120 (tagged "S" for strong field), /3 = 3000 (tagged "M" 
for medium field) and (3q = oo (tagged "W" for weak 
field) . These ISM field models are identified by the codes 
"US", "UM" and "UW" in Table 1. 

In reporti ng their 3D simulation results, Rusz kowski| 
et al. (20071 emphasized that the stabilizing influence of 
their turbulent ICM magnetic field was effective only if 
the field coherence length exceeded the size of a bubble. 
In order to examine the dependence of that result on the 
details of the field model, we also carried out several sim- 
ulations with a simple tangled ICM field model, which we 



identify with the label "T" . While admittedly somewhat 
unrealistic, the model we employed has the advantage 
of having a simple, well-defined coherence length that 
enables us to test directly the importance of this one pa- 
rameter in determining the dynamical role of the field. 
In this model the field lines initially lie in the horizontal, 
y — z plane with 
B x =0 

B y = B (x)f(x, z)[sia(k x ix) sin(k z iz) + cos(k x2 x) cos(k z2 z)] 
B z = B Q (x)f(x,y)[sm(k x ix)sm(k y iy) + cos(k x2 x) cos(k y2 y)}. 

(2) 

In these expressions, the x, y, and z values are mea- 
sured with respect to the grid origin rather than the 



Fig. 2. — Magnetic field configuration used in the uniform ("U") field models. Line brightness is indicative of magnetic field strength. 
The field lines bend to accommodate the bubble inflation region near the middle of the image. Once initialized, a buoyant bubble would 
propagate from the initialization region toward the top of the image. 



bubble origin. Far from the bubble / — » 1, while / 
decreases exponentially towards zero at the boundary 
of the bubble inflation region, so that this ICM field 
once again does not initially penetrate into the bubble. 



The field structure-scale parameters, ki 



i = x,y,z; j 



1,2) were chosen either to make structures large com- 
pared to the bubble inflation region (kij < ir/ri,) (lead- 
ing to a model tag "Tl" for "tangled on large scales"), 
or comparable (fcjj ~ 7r/fb) (leading to a model tag 
"Tb" for "tangled on a bubble scale"). The actual num- 
bers were k x % = 7r/(1.2rb) and k X 2 = n/(0.7rb) for the 
"Tb" models and k x \ — 7r/(6.0?"b) and k x2 = 7r/(3.5n,) 
for the "IV models. The remaining parameters were 
set so that fcjj = (1 + 5)k Xl j ( i=y,z; j = 1,2), where 
S = ±0.04, providing similar scales of field variation. 
In this case we included only a nominally strong field, 
/?o = 8ttP (x) /Bq(x) « 120. Note in this field model that 
the local magnetic pressure varies substantially within 
horizontal planes of constant Pq(x), so local variations in 
(3 are large. To quantify this, we empirically measured 
the average value of (3 in the ICM and found that j3 m 130 
with a standard deviation of comparable magnitude. The 
characteristic field strength once again decreased with 
height. These (strong field) ICM models are identified 
in Table 1 by the labels "T L S" and "T B S". Figure [3] 
graphically illustrates the two tangled field structures. 



2.2. The Bubbles 

While it has been common to conduct bubble sta- 
bility s tudies beginning with static, pre-formed bubble s 
{e.g., Robinson et al. 2004 Ruszkowski et al. 20071, 



instead Jones & De Young (20051 gently "pressure in- 
flated" their bubble structures over a period of at least 
10 Myr during their simulations. Since real relic bubbles 
are not created instantaneously, they argued that rele- 
vant stability issues, including magnetic field geometries, 
may be sensitive to that fact. As a simple compromise 
between static bubbles and full MHD jet simulations to 
create bubbles, they fixed the density and pressure inside 
a localized "inflation region" , allowing their bubbles to 
expand subsonically into the ICM over the inflation time 
period, t,. At the end of that time the inflation region 
was allowed to relax naturally within the surrounding 
flow. We have extended that scheme to 3D for five of the 
computations presented here, and, for direct comparison, 
also carried out seven simulations with full-sized spher- 
ical bubbles pre-formed in the same ICM conditions as 
for the inflated bubbles. 

The bubble inflation region in the simulations reported 
here consisted of a sphere of radius — 2 kpc, cen- 
tered at xt, = 13 kpc, y = z = 0. Conditions in this 
region were kept fixed for times t < U using density, 
Pb = VPo(xb), with i] = 0.01, and twice ambient pressure 
(Pb = 2P (xb))- These conditions were relaxed once t 
reached ti, so that the inflation region was allowed to 
evolve naturally for t > ti. The passive Cf tracer was 
set to unity inside this region when inflation was under- 
way, so that during subsequent evolution gas originating 
with bubble inflation has Cf — 1. We either set U = 
to simulate pre-formed bubbles or ti = 10 Myr for in- 
flated bubbles, which we identify using an added model 
tag "-P . We note that all of the bubble simulations re- 




Fig. 3. — Magnetic field configurations used in the Tg (left) and (right) models. Line brightness is indicative of magnetic field 
strength. Once initialized near the very bottom of the image, a buoyant bubble would propagate from the initialization region toward the 
top of the image. 



ported here involved inflation times U shorter than the 
time necessary for the inflating bubble to expand to the 
size of the local ICM scale height, which is about 25 Myr 
in this case. As |Jones fc De Young| p005[ ) pointed out, 
when bubble inflation is maintained longer than this, the 
rising bubble creates a crude de Laval nozzle in the ICM 
that expels bubble material at mildly supersonic speeds 
upward into the ICM. This results i n the formation of a 
plum e structur e evident in both 2D (Jones & De Young 
|2005[ ) and 3D ( jBriiggen et~aT1|2002| ) simulations. If this 
plume does not develop, then initially spherical bubbles 
comparable in size to the local ICM scale height tend to 
evolve rapidly into a torus around a vertical axis before 
disruptive instabilities can develop. We will discuss this 
behavior in §3.1. 

In most of the simulations reported here the bubble 
inflation region itself was unmagnetized. However, in 
two models, denoted by "-F" , we incorporated a tangled 
internal field derived from the vector potential 



A r (r, 9,4>) — — £?b(Y)rsin 9 sin (mq 
m 



(3) 



defined in spherical coordinates centered at the origin of 
the bubble inflation region with the polar axis along z 
and 4> = along x, the direction of bubble propagation. 
This leads to the circumferential field 

(4) 
(5) 



B = B g 9 - 

with components 

Bg = Bt,(r) sin#cos(m(/ 
and 



Ba, — B(,(r) sin(2f?) sin(r7j<j 

m 



(6) 



where m is an integer. In our simulations m = 2, to 
allow for multiple full field reversals over the azimuthal 
direction. This field has zero magnetic helicity, since 
everywhere A ■ B — 0. We set Bb(r = 0) = and in the 



two "-F" models increased it linearly to Bb(r = 0.8fb) = 
Bo^Xb), whereupon it decreased quadratically to zero at 
r = riy- This maximum value Bo(xb) is approximately 
four times the average ICM field strength and more than 
twice the peak ICM field strength, each measured at x = 
x b - 

3. RESULTS 

We conducted a series of 12 simulations, outlined in 
Table 1, designed to examine the nature of 3D bubble 
propagation and in particular to explore the effects of 
magnetic fields on development of bubble shape and sta- 
bility. Figures [3] -[8] illustrate the structures of each sim- 
ulated bubble at the end of its simulation. These figures 
and their associated animations display volume render- 
ings of the passive color variable C/, to highlight the 
volumes occupied by material that originated within the 
bubble inflation region. To a good approximation this 
traces the regions where the density is substantially be- 
low the neighboring ICM. Figures [3] and [5] show unmag- 
netized bubbles "inflated" into ICMs spanning a range of 
magnetic field strengths and coherence lengths. Figure [6] 
shows the analogous states for "uninflated" bubbles cre- 
ated in the same environments as those shown in Figures 
|4j and [5] Figure [7] displays the forms of bubbles inflated 
into ICMs containing a strong field tangled on different 
scales, while Figure [8] shows two examples of uninflated 
bubbles formed in quasi-uniform and tangled fields. We 
now outline the dominant dynamical developments that 
these images illustrate and the roles of the ambient and 
internal magnetic fields. Comparison of Figures [5] and [6] 
shows that there are only minor differences for a given 
set of other parameters between inflated and uninflated 
bubbles for the short inflation times we used. Primar- 
ily, the inflated bubbles are somewhat larger than the 
uninflated bubbles, as one would naturally expect. The 
basic bubble shapes, howeve r, ar e quite similar. Specif- 
ically, ring structures (see § 3.1l are formed initially in 



all models, reg ardle ss of inflation, and similar bifurcation 
patterns (see § 3.2) are seen in the US and US-I models. 
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Fig. 4. — Volume renderings of the passive "color" variable for inflated bubbles near the end of their evolution, where light colors represent 
bubble material and dark represent the ICM. The ambient magnetic field is along the line of sight. Animations of these quantities as seen 
from several different angles are available at http://www.astro.umn.edu/groups/compastro/ under "Buoyant Bubbles in Galaxy Clusters." 
Top left: UW-I model. Top right: UM-I model. Bottom left: US-1 model. Bottom right: TsST model. 

Given the minor effect of a short inflation period, de- 
tailed discussion of differences between the inflated and 
uninflated bubbles seems unnecessary. 

ft is useful for this discussion to establish up front some 
representative timescales, beginning with r s .o = h/c s .o, 
the characteristic sound crossing time over a scale height 
in the ICM. This is about 15 Myr where our bubbles are 
initialized (x — 13 kpc) and increases approximately lin- 
early with x. The internal bubble sonic timescale would 
be t S) (, = T St o{rb/h)^/rj, much faster than t s . We dis- 
cuss in §3.3 details of the bulk upward motions of the 
bubbles, where we establish their terminal upward ve- 



locity to be typically u t < l/3c Sj o- Thus, the time for 
a bubble to rise across its own radius would be ~ 
r b/ u t = T s,o{ r b/h){c S fi/ut)', this is roughly comparable 
to Ts,o for our bubbles, which initially satisfy r^/h ~ 1/6. 
Timescales for possible disruption from instabilities can 
be estimated from the linear growth rate of perturba- 
tions comparable in size to the bubble radius r^. Using 
the condition r\ <C 1 the hydrodynamical R-T disrup- 
tive timescale would then be t^t ~ \J r b/hr St o, so t^t is 
somewhat less than t s _q or t\, for our bubbles. The analo- 
gous disruption time for the hydrodynamical K-H insta- 
bility, as derived from the dispersion relation for pertur- 
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Fig. 5. — Same as Figure]!] seen from above and with the ambient magnetic field trending diagonally from upper left to lower right. 



bations at a tangential discontinuity (see |Landau fc Lif- 
shitz|19 59 for example) and letting the velocity shear be 
~ u t , would be t kh ~ (r b /h)(c sfi /u t )( l/y/fj)T St0 . Given 
77 <C 1, the K-H disruption time would be significantly 
longer than any of the other listed timcscalcs. 

3.1. Bubble Morphology: Ring Formation 

One striking feature of the early dynamical evolution 
of all the bubbles we simulated is the transition from 
the initial spherical form into a ring or torus around a 
roughly vertical axis. Examination of the velocity fields 
reveals strong flow up through the center of the rings and 
a clear circulation within and around them; i.e., the bub- 
bles evolve into vortex rings or "smoke rings" . Axial sym- 



metry is broken by the ambient magnetic field, although 
the ring-like form is still evident even at late times except 
for the strong magnetic field cases. (Note that "strong" 
in this case still refers to a magnetic field (f3 ~ 100) 
that is initially dynamically unimportant and compara- 
ble to the observationally inferred ICM values) Similar 
evolution is evident in the 3D MHD bubble simulations 



reported by Ruszkowski et al.| fl2007|), the 3D hyd rody- 
namic simulations of Pavlovski et al. (2008 20071, and 
an axisymmetric 2D M HD bubble simulation reported by 
Robinson et al. (2004 1. Although rings were precluded 
i n the 2D p lane symmetric bubble s of |Jones fc De Young| 
( |2005[ ) and |Robinson et al.| ( 2004[ ) , the division there of 
the initial cylindrical bubbles into a pair of counter ro- 



9 




Fig. 6. — Same as Figure [5] except for the models without inflation. Top left: UW model. Top right: UM model. Bottom left: US 
model. Bottom right: TgS model. 



tating line vortices with axes orthogonal to the compu- 
tational plane was analogous. 

Since this behavior seems to be a common element 
of both preformed and briefly inflated buoyant bubbles, 
it is important to identify the cause. So far as we are 
aware, no physical explanation has been offered before, 
although it is straightforward. It is easy to see, in fact, 
that the initial pressure balance between the bubble and 
the ICM leads directly to a strong upward push through 
the bubble center. Much like a smoke ring, the return 
flow around the outside of the bubble sets up the vortical 
flow that leads to a relatively stable ring. The appendix 
outlines two semi- analytic ways to estimate the accel- 



eration of ICM gas into the bubble from below, which 
we have tested quantitatively against ID simulations. In 
short, because the pressure gradient inside a light bubble 
initially close to pressure equilibrium with its environ- 
ment is set by the weight of a denser, external medium, 
the weight of gas inside the bubble itself is not sufficient 
to balance the pressure gradient. The light gas acceler- 
ates upward inside the bubble, creating a rarefaction at 
the lower bubble-ICM boundary. This propagates across 
the bubble in a time t s ,(, <C rj,. 

In response to the lowered pressure inside the bubble 
rarefaction, the lower bubble boundary, or contact dis- 
continuity (CD), accelerates upward at a rate approxi- 
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mately given by (see equations A5 and A7) 



xcd = -g(x)y/p (x)/p b g/s/rj, 



Fig. 7. — Volume renderings of the passive "color" variable for inflated bubbles near the end of their evolution in a tangled ambient field. 
Top: TflS-I model from two different angles. Bottom: T^S-I model from two different angles. Animations of these quantities as seen from 
several different angles are available at http://www.astro.umn.edu/groups/compastro/ under "Buoyant Bubbles in Galaxy Clusters." 

the potential for the bottom of the bubble to catch up to 
the top while the bubble is still close to its initial location, 
(7) which is what we observe to happen in our simulations. 

For the bubble actually to develop into a vortex ring 
the center of the bubble bottom must move upwards 
faster than the bottom regions near the bubble perime- 
ter; that is, the bottom CD must transform from its ini- 
tial concave upward shape into a convex shape. In a hy- 
drostatic ICM that outcome results in equation [7] from 
the fact that \g\/y/rj is generally greatest at the lowest 
point along the CD. The relative difference in xcd grows 
with ri,/h, so vortex rings form more easily when the bub- 
bles are larger. In our simulations, inflated bubbles are 



where g{x) is the local gravity and, again, r\ — Pb/po(xb)- 
Representative timescales for CD motion suggested by 
equation|7]are t C d = y/h/\g\ f]i = (/i/ c *,o) rji = rji t s . , 

and T C D,b = y/ r b/\9\v z = \Zn/h T CD- These measure 
how long it takes the bottom of the bubble to rise up 
through a scale height and across the scale of the bubble, 
respectively. Note that so long as r\ -C 1, TcD,b < T RTi 
and this development within the bubble will be faster 
than disruptive R-T instabilities. If t cD,b < t&, there is 
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Fig. 8. — Volume renderings of the passive "color" variable for bubbles with and without a tangled internal bubble field near the end of 
their evolution. Top left: US model. Top right: US-F model. Bottom left: Tr,S model. Bottom right: Tr,S-F model. Animations of these 
quantities as seen from several different angles are available at |http: / /www.astro. umn.edu/ groups/compastro/ under "Buoyant Bubbles in 
Galaxy Clusters." 



larger compared to uninflated bubbles when they are "re- 
leased" , so they are somewhat more strongly driven into 
the ring transformation. If we ask that the lowest point 
of the bubble CD is driven up to overtake the side of the 
bubble within a time tj, , using equation [7] and assuming 
the acceleration is constant, we can roughly estimate the 
condition for ring formation to be r^/h > r\* (u t /c s fi), 
which is only weakly dependent on r\. For our models, 
this constraint would be r^/h > 1/9. This condition is 
satisfied by our bubbles and als o by the bubbles simu- 
lated by Robinson et al. (2004) and |Ruszkowski et &\\ 
p007|. 



As a concrete illustration of this effect we can con- 
struct a simple 2D planar model assuming a bubble of 
initial uniform density, pg, in an isothermal ICM with 
constant gravity in the —x direction, so that po(x) — 
Poexp(—x/h) and Pq(x) = Poexp(—x/h), with h — 
c s o/(75) — V7- The 7 factor appears because we de- 
fined the scale height, h, using the adiabatic sound speed. 
The gravity in our simulations is not constant, but scales 
approximately as g oc 1/x, as noted above. This actu- 
ally slightly enhances the effects described in this sim- 
pler model. We were unable to find a closed solution for 
the variable gravity case, but verified numerically that 
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the CD acceleration is similar to the constant gravity 
model outlined here. For our simple analytic model we 
parametrize the bubble density to be ps = Tjpo, and de- 
scribe the initial shape of the bottom CD by the curve, 
s(y). It is also convenient to define time scale for each 
point along the CD, 



h r)*y/j Hv) 

t(s) = — e <ifc 

c s .o 

T~CD iM 

= e . 



(8) 



Note that r(s) increases strongly as s(y) increases up- 
wards. Assuming equation [7] applies throughout the CD 
evolution, we obtain an estimate for the position, veloc- 
ity, and acceleration of the CD to time t to be 

"(e* /T 



x(y) = 2/iln 



x(y) = 2 



1) 



21 



at/r 



4e*A 
1 

2*/ r + 1 
P t/r 



(9) 



(e'/T + 1)2 



At the start, when t = and = s(y), the acceler- 
ation is x(y) = (\g\/y/rj)exp(—s(y)/(2h)). By the time 
t = t(s), the point initially at s has risen about half 
the scale height, h. If the bubble is initially spherical 
this timescale is shorter by a factor exp(— n,/(4ft,)) at the 
lowest point on the bubble than at the elevation of its 
center. 

We illustrate this behavior in Figure [9] where the 
evolution described by equation [9] is shown for three 
initially semi-circular, concave upward CDs (s(y) = 

x b — \J r b ~ V 2 )- The three cases vary the bubble ra- 
dius to scale height ratio; namely, left to right h/rb = 2, 
h/r b = 4, and h/r b = 20. The middle case is closest 
to our simulated bubbles. Each case shows multiple CD 
shapes, corresponding to different times in the CD evo- 
lution with a time interval of At ps 0.08tcd between 
consecutive shapes. Except for the third case, where the 
scale height is much larger than the size of the bubble, 
the CD is quickly inverted to a convex upward shape, 
consistent with what is seen in the simulated bubbles. 
While this model is a simplification of the real dynamics 
in 2D or 3D, it does demonstrate the basic initial physical 
behavior of the bottom of the bubbles. 

As a result of this behavior the ambient gas drives up- 
wards in the center of the bubble, while the buoyant 
rising of the bubble body establishes a downward flow 
around the outside, thus generating a vortex ring. Our 
previous comparison with the growth rate for disruptive 
R-T instabilities showed that this ring development is 
faster when the bubble density contrast is large. We 
expect, then, any light, initially quasi-spherical bubble 
with h < rb to be dominated by this effect, so long as 
the bubble inflation process is terminated before the bub- 
ble top has risen through a scale height, or roughly an 
interval several times r S: o = h/c S fi. As mentioned previ- 
ously, when the inflation process is maintained over inter- 
vals longer than this a de Laval nozzle forms and bubble 
material expands into the upper atmosph ere at speeds 
~ c s ^b ~ c S fl/^fq ( Jones fc De Youiig||2005 1 . This flow is 



too fast to be overtaken by the lower boundary moving 
upwards at a speed slower by a factor ~ 771 according to 
equation [9 If the bubble inflation takes place at the cen- 
ter of the ICM, if the initial bubble was large compared 
to the local scale height, or if it involved a driven flow by 
a jet from outside the bubble, the n this dynamical be 
havior would not be releva nt (e.g. 
IBasson & Alexander||2003| lOmma et a 



Binney|2004j |O'Neillirjones|2009|). On the other hand 



Reynolds et al. 2002 

' ||2oo4| |OmmT" 



Favlovski eFaXl j 2008 2007|) have n oted that radially- 
directed Kutta-Zhukovsky forces (see Landau fc Lifshitz] 



19591 f° r example, for a discussion of Zhukovsky's theo- 
rem) may characteristically transform flattened bubbles 
into vortex rings regardless of the details of inflation. 
This is distinct from the process we have identified, how- 
ever, in that Kutta-Zhukovsky forces require the presence 
of both circulation and bubble asymmetry to produce a 
vortex ring. In contrast, the mechanism of ring formation 
that we describe operates immediately, even on initially 
spherical and motionless bubbles, so long as their pres- 
sure profiles are set by the equilibrium structure of the 
comparatively dense ICM. 

We close this portion of the discussion with two simple 
observations about rings formed in unmagnetized media. 
The first is that, in our models the rings in the UW and 
UW-I simulations had roughly steady tube radii, with 
r ~ r b . The tube radius was larger in the UW-I case, 
basically because it was bigger when "released" . Sec- 
ond, in the absence of magnetic fields the ring-like bubble 
structures formed in our simulations were relatively sta- 
ble over the simulation times, although they developed 
some bending mode oscillations, probably from the Wid- 
nall i nstability flWidnall fc Sullivan|[l973{ |Widnall et al 



1974) ) that is known to affect vortex rings, and surface ir 
regularities likely due to K-H instabilities. The observed 
coherent behavior is expected for vortex rings when the 
effective Reynolds number, R e is less than about 10 3 , al- 
though for much higher Reynolds number flows vortex 
rings are seen ex perimentally to be unstable to turbulent 
break-up (e.g., Maxworthy 1972). It is not possible to 
define a rigorous Reynolds number for calculations with 
our inviscid MHD code; however, based on code tests 
and the fact that the tubes were spanned in our simu- 
lations by ~ 30 cells, we can estimate an effective (fast 
mode) Reynolds number for flow around the ring to be 
R e ~ 500 - 1000 {[Hyii_eial.||20()()|. The 3D simulations 



by Ruszkowski et aT. 



with a Reynolds numbe 



POTT 

ber, R, 



applied an explicit viscosity, 



2000. The ring structure 
seen in their unmagnetized bubble simulation appears to 
have survived intact for w 100 Myr, although those au- 
thors also note the presence of small scale K-H instabili- 
ties. Thus, if such rings are indeed destructively unstable 
on timescales of interest, substantially higher resolution 
than any currently available seems to be necessary to 
capture that behavior. 

3.2. Bubble Morphology: Magnetic Influences 

The early evolution of the bubbles (i.e., ring forma- 
tion) was essentially independent of the magnetic field 
properties we applied. That is expected, since the initial 
Maxwell stresses from even the strongest ICM field, with 
(3 = 120, were insignificant. That result mirrors our 2D 
simulations as well. We note, however, that |Ruszkowski| 
et al. (2007) found varying degrees of ring formation be- 
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Fig. 9. — Semi-analytic model for the evolution of the lower surface of a buoyant bubble for varying ICM scale heights. Left: h = 2r^. 
Center: h = 4r&. Right: h = 20r^. The two axes are 2D spatial dimensions in units of the bubble diameter, where +x is the direction of 
bubble propagation. Time increases from bottom to top, with a time interval of At O.OStcd between consecutive CD shapes. 



tween their magnetized and non-magnetic models. While 
it is clear that all of their models, like ours, are sub- 
ject to the initial forces that would cause ring formation, 
their magnetized models featured initial fields ((3 ~ 40), 
both internal and external, that were apparently suffi- 
ciently strong to prevent the completion of this process. 
Although we have not attempted here to model fields 
strong enough to disrupt ring formation, our simulations 
do show that the inclusion of even a moderate ambient 
magnetic field ((3 = 3000) can have a clear impact on 
evolution of the ring structures, provided the coherence 
scale for the field exceeds the size of the bubble. 
Consider first the simpler cases with quasi-uniform am- 



bient fields. The "strong" field cases, US, US-I, with 
= 120 show this behavior most obviously, so we exam- 
ine those first. The lower-left images in Figures [4] and 
[H and their associated animations illustrate that these 
bubbles developed a wishbone shape; that is, the bubbles 
were bifurcated. The initial ring became elongated in the 
horizontal direction orthogonal to the ICM field. As the 
animations show, the upward motion of the central por- 
tion of the ring aligned with the ICM field was stalled and 
then reversed, while two plumes of bubble material on 
opposite sides of this bifurcation continued to rise. This 
behavior results when the external field gets trapped in 
R-T or the aforementioned Widnall instability-induced 
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Fig. 10. — Field lines in the US-I model near the end of the simulation time, showing the magnetic field structures that have bifurcated 
the bubble. The color on the lines represents field strength (lighter=stronger). 



oscillations of the upper boundary of the bubble. Ini- 
tial upward motions stretch and amplify those field lines, 
which are anchored in the ICM away from the bubble, 
until the magnetic tension is sufficient to resist contin- 
ued upward motions. Eventually, this magnetic tension 
exerts a force that exceeds the buoyant force, causing 
those portions of the bubble in which the ICM field is 
embedded to reverse direction and sink. R-T and K-H 
instabilities in the plane of the field and the bubble mo- 
tion (the x — y plane) become inhibited, of course, even 
though perturbations orthogonal to the field are not. We 
verified that in these regions can approach unity and, 
more importantly, that the Maxwell stresses (i.e., j x B) 
are comparable to and sometimes dominate the gas pres- 
sure and gravitational forces (VP and pg) there during 
this period. This highlights the fact that the (3 parame- 
ter is not sufficient to evaluate the role of the magnetic 
field, since the relative stresses depend not only on the 
magnitude of the magnetic field, but also on the scale 
and direction of its variation. On the other hand, since 
bubble motions in the x — z plane are not magnetically 
inhibited, portions of the bubbles slide through the field 
lines, creating the plumes. The magnetic field lines and 
intensities during this stage of the evolution of model 
US-I are illustrated in Figures 10 and 11 respectively. 



This evolution is qualitatively similar to that seen for 
dense supers onic, spherical clumps modeled by |Gregori| 



et al. (1999), where an initially uniform ISM magnetic 
field wrapped around the clump, eventually also bifur- 
cating it. Like our 3D bubbles, the boundaries of the 



Gregori et al. (1999) clumps w ere relatively s harp. On 
the other hand, the outcome of |Gregori et al. ( 1999 1 ap- 



pears differ ent in some ways from the simulation results 
reported by Dursi & Pfrommer (2008) of dense, centrally 
condensed, spherical clumps moving from an unmagne- 
tized environment into a medium containing a laminar 
"draping" magnetic field. The latter computations were 
designed to study formation of cold fronts during cluster 
mergers. Those clumps developed a strongly magnetized 
layer on their noses, but at least on timescales required 
for them to sweep by their own masses they were not 
bifurcated. 

Bifurcation behavior seen in the US and US-I simula- 
tions is also visible in the UM and UM-I models, although 
it is much less pronounced, due to the much weaker ini- 
tial fields. In the upper-right images in Figures [4] and [5] 
and their associated animations, we can see that vortex 
shedding leads to a buildup of material in the direction 
orthogonal to the uniform field. This asymmetry appears 
to be caused by the same variety of magnetic stresses 
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Fig. 11. — Midplane slices showing the ratio of the magnetic field strength to its initial value (B/Bq) for the US-I model in the x — y 
(left column) and x — z (right column) planes. The top row shows t = 25 Myr, middle t = 50 Myr, bottom t = 100 Myr. The colorbar 
ranges from [0:6]. 
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that are present in the US and US-I simulations, and 
over longer times it is likely that these rings would also 
become bifurcated. It is also clear that the portion of the 
UM bubble oriented along the field lines is less elevated 
than the rest of the structure, which is also similar to the 
behavior of the bifurcated bubbles. 



As emphasized by Ruszkowski et al. (20071 the influ- 
ence of the magnetic field on bubble evolution depends 
also on the geometry of the field. In particular, if the 
ambient magnetic field is tangled on scales less than or 
comparable to the scale of the bubble, then field lines 
are not effectively stretched and amplified as the bubble 
lifts upwards, in comparison to what happens if the field 
structure is laminar. We explored that issue further in 
our experiments using our simpler tangled strong field 
models, T^S and T^S outlined in §2.1. In short, these 
simulations have ICM field strengths comparable to the 
US models, but field lines that are tangled on scales of the 
initial bubble circumference and a factor of five greater, 
respectively. Figure [7] shows a comparison of the bub- 
ble morphologies for the T^S-I and T^S-I models, while 
comparisons with uniform field cases are available in Fig- 
ures gj [5] |6j and § 

In both the T^S-I and T^S-I cases the ring structures 
that form remain intact to the end of the simulation, 
whereas for the analogous US-I simulation the ring was 
bifurcated, as discussed previously. For the smaller scale 
field structure (T^S-I), in fact, the bubble ring is hardly 
different from the unmagnetized case. An examination of 
the magnetic field properties explains this simply, since 
it turns out the field is generally not amplified compared 
to initial conditions. In fact, as a result of reconnection 
that takes place in the flow around the bubble, the (3 of 
the ICM plasma near the bubble is actually increased. 

The situation is more complex in the T^S-I simula- 
tion. While the ring structure survives to the end of the 
calculation, there are localized regions where the field is 
strongly amplified by stretching in response to the bub- 
ble's rise. This leads to substantial distortions in the 
bubble, but in the absence of coherence in the field over 
very large scales the dynamical impact of the Maxwell 
stresses is irregular. It is likely, however, that over still 
longer timescales this bubble would be disrupted by the 
magnetic field. Our results on the effects of field tan- 
gl ing are, therefore, i n q ualitative agreement with those 
of iRuszkowski et al.| ( [2007| . 

Having noted that magnetic reconnection can increase 
the value of (3 by reducing field strengths near the bub- 
ble in our models of tangled ICM fields, it is worth 
briefly detailing the nature of this process in our sim- 
ulations. At the most basic level reconnection is a topo- 
logical transformation of the magnetic field made pos- 
sible by local d issipation where adjacent fi elds reverse 
or braid {e.g., Lazarian & Vishniac 20081, setting up 
concentrated currents. Like other general purpose MHD 
codes our numerical algorithms include no explicit pre- 
scription for reconnection, per se; rather, reconnection 
is enabled by the effects of numerical dissipation where 
fields reverse or braid close to the resolution scale of the 
computational grid. Although numerical in nature, this 
process includes essential physics involved in reconnec- 
tion; namely, it exactly conserves energy and momentum, 
and maintains a divergence free field, while allowing such 
physical, reconnection-related, nonideal MHD processes 



as tearing modes to develop. Consequently, the code in- 
cludes the effects of physical reconnection by canceling 
field and channeling field energy and stresses into ther- 
mal and kinetic energy and momentum appropriate to 
the details of the flow and the field. In our models of 
uniform ICM fields, we expect reconnection to play no 
significant role since localized field reversals are rare (see 
the field configuration in Figure 10 for example). In 



our models of tangled ICM fields, however, reconnection 
surely does take place, since local reversals and braiding 
are common. Reconnection, therefore, modifies the fields 
as the system evolves. While this changes the detailed 
structures of the fields, it does not substantially modify 
field structures on scales beyond the coherence length of 
the field. 

Finally, we touch on the properties of the two simula- 
tions with magnetic fields internal to the bubbles (cases 
US-F and T L S-F). As introduced in §2.2, those calcula- 
tions were initialized with a tangled magnetic field con- 
tained within the bubble. The maximum magnetic field 
just inside the bubble boundary was approximately a fac- 
tor of four greater than the average external field at a 
comparable ICM height. As it turns out, the internal 
field quickly dissipated as the bubbles expanded due in 
part to weakening of the field from flux freezing, but 
also in response to reconnection coming from the tan- 
gled nature of the original field geometry. Consequently, 
the presence of internal magnetic fields was almost com- 
pletely irrelevant to the evolution of our models. Fig- 
ure [8] provides a comparison between the US and T^S 
(left) and the US-F and T L S-F (right) models. There 
are subtle differences present in each pair, but the inter- 
nal fields clearly have not played significant roles. We 
found a somewhat different conseq uence in our 2D sim- 



ulations ( Jones fc De Young] 2005[ ) . In those simulations 
circulations established inside the bubbles stretched and 
amplified the internal fields, which, however, were ini- 
tially circumferential rather than tangled. This points to 
the importance of the field geometry and to differences 
between 2D and 3D evolution. In 3D bubbles that were 
initially turbulent, so that turbulent field amplification 
was ongoing, we might once again expect some enhance- 
ment of the field so long as the turbulent motions were 
strong enough to control the bubble field evolution on 
timescales of order r Si o. 

3.3. Bubble Upward Propagation 

Besides morphology the most applicable measure of 
bubble dynamics is a bubble's upward motion within the 
ICM. Figure [T2| illustrates this behavior for the bubbles 
evolving in the quasi-uniform external fields. The quan- 
tity plotted is the mass-weighted mean elevation of the 
bubbles, x av = J CfxpdV/ f C/pdV. As the discussion 
of the previous subsection suggests, the upward motions 
of the bubbles placed in a tangled external field are simi- 
lar. The motions indicated in Figure[l2]are close to those 
found for the analo gous 2D bubbles simulated by [Jones 
& De Young ( |2005[ ). They indicate an upward velocity 
that is initially about 0.3c Sj0 , but which slowly decreases 
over time, as the bubbles rise higher in the ICM. The 
inflated bubbles decelerate somewhat less than the unin- 
flated bubbles. 

The upward bubble motion is, of course, driven by 
buoyancy and limited by dynamical drag and Maxwell 
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Positions of the mean heights of bubbles in uniform ambient magnetic fields. Top: Inflation models. Bottom: Models without 



stresses. Along the way local dynamical features, such 
as R-T and K-H instabilities, as well as buoyancy-related 
dynamical evolution of the bubble-ICM interface influ- 
ence the sh ape and integrity of the bub bles, as we have 
discussed. |Jones fc De Young (20051 found that the 
upward terminal velocity of their 2D bubbles were, in 
fact, well-described by the simple balance between the 
buoyant force, Ft, uoy w |si/>oH 5 and the drag force, 
Fd ~ CdAbPoul, where Vb and Ab are the bubble vol- 
ume and "cross sectional' area, respectively, Ub is the 
bulk speed of the bubble, Cd ~ 1 is the drag coefficient, 
and we utilized the condition r) = pb/po 1. Initially, 
when the bubble is quasi-spherical Vb ~ (4/3)7rr^ and 



Ab ~ 7rr^, while after ring formation Vb ~ 2ir Rr^ and 
Ab ~ AirRrb, where R is the radius of the ring. Balancing 
the buoyant and drag forces gives the familiar expression 
for the terminal upward velocity, 



( \g\Vb 

\C d A b 



(10) 



The factor Vb/Ab is (4/3)rb for a sphere and (n/2)rb for a 
ring, so in either case the terminal velocity can be written 
to sufficient accuracy simply as 



ut 



(11) 
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where we set \/Vb/(AbC'd) = 1 and utilized the definition 
h = <? sQ /\g\. The time to reach this equilibrium speed is 
very brief, about 



n 



rjh 



Jones & De Young ( 2005 ) computed Vf, and At numer- 
ically from the simulation data, finding a close match 
between the actual upward velocity and that predicted 
by equation |11| For our purposes here a simpler anal- 
ysis suffices. Applying the initial bubble conditions to 
equation [TTj gives ut ~ 0.4c s ,o, which is reasonably close 
to the measured velocities. The slow deceleration that 
takes place at higher elevations can be understood from 
the fact that the tube radius in the ring r^, does not 
increase substantially as the bubble rises. The radius 
of the torus, R, does increase as the bubble rises, but 
that has little impact on the terminal velocity, which is 
proportional to yj r^/h. Since h increases in our simula- 
tions with elevation, the terminal velocity is smaller at 
higher elevations. The inflated bubbles are less deceler- 
ated than the uninflated bubbles, because their tube radii 
are larger, as mentioned before and as can be seen in the 
figures. In the first approximation the upward trajecto- 
ries in Figure 12 are independent of the ambient magnetic 
field. There is, however, a notable dip in the curves for 
the cases US, US-I and US-F beginning around t ~ 30 
Myr. This corresponds, of course, to the stage during 
which the ambient field bifurcates those bubbles. Some 
of the bubble matter trails behind and actually begins to 
fall, as a result of the reaction of the stretched field lines. 
At the same time, however, bubble plumes continue to 
rise. The mass- weighted bubble velocity decreases ac- 
cordingly, but eventually recovers as the plumes, which 
have characteristic sizes greater than the ring tube size 
at the time of bifurcation, move upward more rapidly 
than the unbifurcated rings. 



4. CONCLUSIONS 

We conducted an ensemble of simulations designed to 
explore the 3D nature of buoyant bubbles propagating in 
the magnetized ICM. 
Our most important results are summarized here: 
1. Consistent with previous 2D and 3D simulations 
we confirm that magnetic fields of strength (3 ~ 120 can 
substantially affect the morphology of buoyantly rising 
bubbles. The influences depend significantly on the ge- 
ometry of the field. In general, the results here do not 
show the extreme bubble break up and subsequent mix- 
in g that are seen in th e purel y hy drodynamic simulation s 
of |Bruggen fc Kaiser| ( |2001| ) and [Briiggen et al.| ( |2002| ), 
which is in accord with earlier MHD calculations! In 
uniform horizontal fields, rising bubbles can become bi- 
furcated as a result of magnetic field trapping along the 
upper contact discontinuity and subsequent stretching 
across the center of the bubble. If the ambient field is 
tangled on scales smaller than the size of the bubble, 
the dynamical influence of the field is greatly reduced, 
since those fields are not effectively amplified by the up- 
ward motions of the bubbles. If the field tangling scale 
is significantly larger than the bubble, the dynamical in- 
fluence of the field can be significant, although relative 
to a uniform field of comparable strength the influence 



is retarded and less coherent than that due to a uniform 
external field. 

2. If the ambient fields have coherence scales larger 
than the bubbles, the upward bubble motions can am- 
plify the field strengths locally and transport magnetic 
flux to higher elevations in the ICM. This field transport 
and amplification may lead to locally enhanced mixing 
and lifting of the ambient ICM, and this process can in 
principle contribute to the suppression of cooling flows. 
However, the localized nature of this process does not by 
itself provide an effective means of raising the ICM tem- 
perature in the core on the global scales that are needed. 

3. The buoyant motions of bubbles can be approxi- 
mated to zeroth order by a simple balance between buoy- 
ant and drag forces, but the presence of magnetic fields 
and their consequent dynamical influence do have a clear 
effect upon the upward motion of the bubbles, as is de- 
scribed in Section 3.3. 

4. Hydrostatically formed bubbles, and in particular 
bubbles inflated slowly in the ICM, provide a useful tool 
for understanding the interactions of buoyant material 
with ambient magnetic fields. In particular such models 
are applicable to the evolution of the relic radio bubbles 
seen in the ICM of many rich clusters. In general, such 
spherical bubbles having scales comparable to the scale 
height of the local ICM tend to produce toroidal bubble 
morphologies. If, on the other hand, the objective is to 
understand the morphology of radio bubbles directly in- 
flated at the termination point of AGN jets, then it is 
necessary to evolve that structure from the dynamics of 
a jet rather than to artificially create a bubble in hydro- 
static equilibrium with the ICM. 
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APPENDIX 

ACCELERATION OF THE BOTTOM BUBBLE-ICM INTERFACE 

We outline two simple ways to estimate the early ID acceleration of the lower contact discontinuity, CD, between 
a buoyant bubble and its environment when the bubble is formed in an initially hydrostatic medium. Gravity points 
downwards in the — x direction. The motion of the CD is a response to a strong rarefaction that forms inside the 
bubble above the discontinuity. The rarefaction takes the ambient gas below the bubble out of hydrostatic equilibrium, 
so that this gas accelerates upwards towards the top of the bubble. 

The initial pressure distribution everywhere is assumed to correspond to hydrostatic equilibrium in the ambient 
medium; namely, 

dPa(x) , . . . 
d = Po{x)g{x) < 0, 

where g(x) is the signed gravitational acceleration. The initial CD is at x 
p(x > s) = pb <C Po = p(x < s). As in the text we define r\ = pb/po <C 1. 
The net acceleration of gas is 

p ox 



(Al) 

s, with the bubble on top, so that 



(A2) 

From equation |Al| 'x = initially for ambient gas, but not for the much lighter bubble gas. The net acceleration of 
the bubble gas is, instead 



Po 
Pb 



1 



V 



(A3) 



The acceleration of the underlying ambient material occurs in response to a rarefaction that forms between the two 
substances, increasing the local pressure gradient above the initial hydrostatic value. The rarefaction propagates at 
the local sound speed, which is c s ^ in the bubble and c Si o in the ambient medium. Because of the initial pressure 
equilibrium, c s _q = y/rj c Sj b *C c s ,b- After a short time, At, the rarefaction moves into the bubble a distance, Aa; r ,t> = 
c s ^At, and oppositely into the ambient medium a smaller distance, Ax r ^ = — c S fiAt. The pressure inside the 
rarefaction, P r , is approximately the initial pressure just above the rarefaction, so P r « P(x — s + Ax r b)- Using 
this gives P r rs P(s) + pa gAx r b- The pressure difference across the rarefa ctio n propagating into the 
:ium then becomes approximately AP sa — po g(Ax r> b — A:r r) o) Using equation A2 the net acceleration of 



equation 



Al 



ambient medi 
ambient gas in the rarefaction is then approximately 
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Po 
9(x) 



-p g(x)(Ax rt b " Ax r fi) 



lAavnl 



Ax r i 



Ax r o 



(A4) 



The ratio Ax r> b/ Ax r fl is simply the negatively signed ratio of the sound speeds, — c Sj b/c s ,o, so associating this motion 
with the CD gives 



xcd 



-sO) 



c s ,b 

C s ,0 



?7 2 



(A5) 



Alternatively, one can estimate the motion of the CD from the solution of an approximate Riemann prob- 
lem representing the initial conditions. Ignoring gravity for now and setting 7 — 5/3, the Riemann invariant, 
J_ = u — 3c Si b = constant, applies inside the rarefaction propagating into the bubble from its lower boundary, 
while inside the associated rarefaction moving into the ambient medium J + — u + 3c s = constant applies. The local 
velocity is given by u. Just above the rarefaction penetrating the bubble the velocity is, according to equation |A3| 
Ub ~ — gAt(\ — rf)/r), whereas just below the rarefaction penetrating into the ambient medium uq = 0. The Riemann 
invariants allow us to find the matching conditions at the CD, and specifically its velocity, 



ucd = u + c s , B 



Uq - u b 



_C S ,0 + c s,b 

Since 77 <C 1 , we have once again (letting 'xcd ~ ucd I At) 



gAt 1-r] 

7J2 1+7^2 



XCD 



(A6) 



(A7) 



We have verified that this behavior is consistent with ID simulations of a contact discontinuity set up with these initial 
conditions. 
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